Critical scaling of jammed system after quench of temperature 
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Critical behavior of soft repulsive particles after quench of temperature near the jamming tran- 
sition is numerically investigated. It is found that the plateau of the mean square displacement of 
tracer particles and the pressure satisfy critical scaling laws. The critical density for the jamming 
transition depends on the protocol to prepare the system, while the values of the critical exponents 
which are consistent with the prediction of a phenomenology are independent of the protocol. 
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I. INTRODUCTION 

The jamming transition has attracted many physicists 
since Liu and Nagel indicated its similarity to the glass 
transition In naive sense, the glass transition is char- 
acterized by a divergence of time scale on a temperature- 
density plane 0-13|; while the jamming transition is an 
athermal transition on a density-load plane as the emer- 
gence of rigidity for materials such as granular materials, 
forms, and colloidal suspensions [5|-|7|. For frictionless 
particles, it is known that the pressure, the elastic mod- 
ulus, and the characteristic frequency for the soft mode 
continuously emerge above a jamming transition point, 
while the coordination number changes discontinuously 
at the point [sl-fTTj. 

The jamming transition was discussed mainly on the 
axis of the density without any load and temperature 
in some pioneer works [sl-fllj. It is instructive, however, 
that critical properties have been clarified when we look 
at the behavior of the jamming on the density-load plane 
in the zero load limit. For instance, critical scaling laws 
for the rheological transition, similar to those in contin- 
uous phase transitions, have been observed for sheared 
frictionless systems [l^ - |2l| , while the discontinuous tran- 
sition and the hysteresis loop are observed in the pressure 
and the shear stress for sheared frictional granular mate- 
rials [H,!!!. 

Coming back to the original idea in Ref. we 
can also discuss the jamming transition on the density- 
temperature plane, i.e. without any load. This approach 
has an advantage to clarify the relationship between the 
glass transition 0-13 and the jamming transition 0-01 j 
because the glass transition is originally defined only on 
the temperature-density plane. So far the behavior on 
the density-temperature plane in the zero temperature 
limit has been studied by some researchers, but the sit- 
uation is still confusing. Indeed, some indicated that 
the critical fraction (j)Q for the divergence of the relax- 
ation time in the zero temperature limit is identical to 
the critical point 0j for the jamming transition of the 
athermal materials 0, , while the others suggested 
that two transition points are different [26l - l30j . It is re- 



markable that Berthier and Witten have numerically con- 
firmed from their simulation for soft repulsive particles 
at low temperature that (i) the relaxation time around 
the glass transition point (j}Q satisfies a scaling relation, 
and (ii) 0g is lower than that for the jamming point 
[29l . [30| . It is also noticed that the separation between 
4>G and 4>} is clearly demonstrated from a simulation for 
sheared soft spheres in the zero temperature and zero 
shear limits [31| . 

Recently, the critical behavior of repulsive particles on 
the density-temperature plane near the jamming transi- 
tion at zero temperature and high density has been stud- 
ied both numerically and theoretically |32h35| |. It is no- 
table that the replica theory gives a prediction on both 
the critical fraction and the critical exponents [s^, 'ss'] . 
The validity of their prediction for the critical behavior 
of the pressure, the energy, and the divergence of the 
first peak of the radial distribution function have numer- 
ically verified JB2, 33]. However, it is unclear whether the 
critical exponents are unique because they might depend 
on the protocol to prepare the system as for the critical 
density of the jamming transition (36l[37|. 

In this paper, to clarify critical behavior on the density- 
temperature plane in the vicinity of the jamming transi- 
tion point, we numerically investigate the value of plateau 
(VP) of the mean square displacement (MSD) of tracer 
particles and the pressure of soft repulsive particles after 
quench of temperature and demonstrate that the criti- 
cal exponents for critical scaling laws does not depend 
on the protocol, while the protocol dependence exists in 
the critical fraction. We should note that the pressure 
for soft spheres [l^l and MSD for hard spheres [S] have 
been numerically measured in the previous papers, but 
this paper is the first report on the numerical study of 
MSD for soft spheres. 

The organization of this paper is as follows. In the next 
section, we will explain our set up and models. In Sec. 
mil we show the results of our simulation on the critical 
behavior for VP and the pressure of the quenched soft 
particles. In Sec. llVi we will show the jamming tran- 
sition density depends on the protocol to prepare the 
system. In Sec. |Vl we will present scaling laws for the 
plateau and the pressure, and theoretically determine the 
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critical exponents. In Sec. I VII we will discuss and con- 
clude our results. 



II. SETUP AND MODEL 

We study a three dimensional system consists of N 
soft spherical particles with mass m enclosed in a peri- 
odic cube of linear size L. Note that the box size L is 
fixed for the most cases, but is changed when we will 
determine the jamming point in Sec. IV. We prevent 
the system from crystallization by using a 50:50 binary 
mixture of spheres of diameter ratio 1.4 which is numer- 
ically confirmed from the radial distribution function, 
where sharp peaks characterizing crystallization do not 
exist d, H, H^, [S^l ■ It should be noted that the critical 
behavior for jamming transition of granular particles is 
unchanged even for a mono-disperse system or a poly- 
disperse system with equal number of particles of diam- 
eters (To, 0.9(70, 0.8(To, and O.Tctq [sl-Hflst. 

For later convenience, let us use dimensionless quan- 
ti ties sca led by ctq for the length, m for the mass, and 
a/ maQ/e for the time, respectively, where we have intro- 
duced a characteristic energy scale e. We assume that 
the interaction between i and j particles is described by 
a pair wise potential 



V{r,j) = (1 - nj/(T,jfe{a,j 



(1) 



where 6{x) is 
e{x) = 1 for 



the Heaviside step function satisfying 
X > Q and 0{x) = for otherwise, 
ij ~ 



i^ij — ~ ''jl and aij = (ai + cTj)/2 with the position r 



and the diameter at of the particle i. 

We start from an equilibrium state at an initial tem- 
perature Ti and a volume fraction (j). Then, we quench 
the system directly to a final temperature Tp, and the 
system subsequently evolves at Tp by the velocity rescal- 
ing thermostat. We use the system size N = 1000. We 
have checked the critical exponents do not change when 
we use N = 4000. We adopt the leap-frog algorithm 
with the time interval At — 0.01. We have verified that 
the choice of the algorithm does not affect the average 
values of the pressure and MSD within the numerical ac- 
curacy when we use the velocity Verlet algorithm and 
At = 0.001. We believe that the initial state is suffi- 
ciently equilibrated. Indeed, as long as we have checked, 
MSD exceeds 10 and we could not find any aging effects 
during the equilibration process. This system has been 
well studied in the previous papers on the ener gy, the 
pressure, and the radial distribution function (32l l33j| . 



III. MEAN SQUARE DISPLACEMENT AND 
PRESSURE 



In this section, we summarize the results of our simu- 
lation on MSD and the pressure. Note that the system 
has a fixed volumed fraction (p or a, fixed volume in this 
section. 
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FIG. 1: MSD (r2(t)) as a function of time t for <?!> = 0.7, 
Ti = IQ-^ and Tf = 10"^ with tu, = 10^ 10*, and 10^ 
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First, let us consider the mean square displacement of 
larger tracer particles at Tp 



NI2 



(2) 



where rL,i(t+tu;) and tw are the position of the larger par- 
ticle i and the waiting time, i.e., the time elapsed after the 
quench, respectively. Here, we ignore the displacement 
of the smaller particles in order to eliminate the effect of 
rattlers @ . The bracket denotes an equilibrium ensemble 
average over the initial configurations. In Fig. [1] we plot 
MSD {r'^{t)) as a function of the time t with (j) = 0.7, 
Ti = 10-2, Tp = 10"^ for the waiting time = 10^, 10'', 
and 10^, where MSD exhibits clear plateaus. The time to 
escape from the plateau increases as the waiting time 
increases, which indicates that the system does not reach 
an equilibrium state within the time window explored in 
our simulation [s^. However, we should note that VP is 
independent of the waiting time tw- 
in Fig. [21 we plot MSD as a function of the time 
t divided by the "thermal" time tt = l/^/Tp for (j) — 
0.62, Ti = 10-2, tw = 10^ with Tp = lO-'', 10-^ IQ-^, 
and 10-^. Thanks to the introduction of the scaled time 
t/r-r, MSD for (j) = 0.62 converges to a master curve, 
which indicates that VP is almost independent of the 
final temperature Tp. For relatively low density case, it 
is known that the particles behave as a hard sphere liquid, 
in which the dynamics is independent of the temperature 
if the time is scaled by the thermal time [1^ [13] . This is 
the reason for the scaling behavior as shown in Fig. [2l 

On the contrary, MSD strongly depends on Tp for 
denser cases. In Fig. [31 we show {r'^{t)') scaled by Tp 
as a function of the time t for <j) — 0.70, Tj — 10~^, and 
tw = 10^ with Tp = 10-4, 10-5, and 10-<^. MSD {r^{t)) 
scaled by Tp converges to a master curve, which indicates 
that VP is proportional to Tp. For this particle is 

completely trapped within a cage and fluctuates around 
its equilibrium position. The reason why VP is propor- 
tional to Tp can be understood as follows. The energy SE 
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FIG. 2: MSD {r^{t)) as a function of the time t scaled by 
the thermal time tt for = 0.62, Ti = 10"^ and t„ = 10^ 
with Tp = IQ-'', 10-^ 10-^ and 10"^ 



FIG. 4: (Color online) The value of plateau rUp as a func- 
tion of the Tp for Ti = 2.0 x 10"^ and = 10^ with 
(j> = 0.62, 0.64, 0.65, 0.66, 0.68, and 0.70. 
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FIG. 3: MSD {r^{t)) scaled by Tp as a function of t for 
(j> = 0.70, Ti = 10~^ and U = 10^ with Tp = 10"*, 10~^ and 
10"^ 



FIG. 5: (Color online) The pressure p as a function of the 
final temperature Tp for Ti = 2.0 x 10"^ and t» = 10^ with 
= 0.62, 0.64, 0.65, 0.66, 0.68, and 0.70. 



due to the fluctuation of its position Sr may be approxi- 
mated as 6E (X \6r\'^. If we assume that the distribution 
of 6r satisfies p{6r) cx exp{—SE /Tp), p{Sr) depends on 
through \Sr\'^/T-F. If we also assume that VP is scaled 
by the size of the fluctuation (|(5rp), it is reasonable to 
obtain the scaling relation as in Fig. |3l 

Here, let us introduce nip as {r^{t)) &t t = tt- We 
should note that {r^{t)) changes less than 10 % for t > 
Tt ■ Figure [5] exhibits rup as a function of Tp for = 
0.62, 0.64, 0.65, 0.66, 0.68 and 0.70 with Ti = 2.0 x lO^^. 
As we have noted, nip is a constant for lower densities and 
is proportional to Tp for higher densities. It is notable 
that rUp behaves as a power-law function of Tp around 
(I) = 0.65 [1]. 

The similar critical behavior can be observed for the 
pressure at the final temperature Tp 

where is the momentum of the particle i and f{rij) = 
— VlAvij) is the potential force. Figure [S] shows the pres- 



sure p as a function of Tp for Tj = 2.0 x 10^^. For 
(j) = 0.62, p is almost proportional to Tp which is one 
of characteristic behavior of hard sphere liquids. On the 
other hand, p is a constant at higher volume fraction 
such as (j) = 0.70, because the pressure is determined by 
the rigidity of contact network of particles. It is reason- 
able that the rigidity of the network is insensitive to the 
temperature near T = 0. 

IV. PROTOCOL DEPENDENT CRITICAL 
FRACTION 

In this section, let us determine a critical fraction 0j 
after the quench. It should be noted that is deter- 
mined not by a simulation under a fixed volume but by a 
simulation by a compress or an expansion of the volume. 

First, we prepare an equilibrium state of the volume 
fraction and the initial temperature Tj. Second, we 
quench the system directly to Tp and keep the temper- 
ature by the velocity rescaling thermostat. Third, we 
further relax the system to the nearest potential energy 
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FIG. 6: (Color online) The critical fraction <^,j((/), Tp, Ti) as a 
function of (j) and Tp for the initial temperature Ti = 2.0 x 
10"^ 




FIG. 7: (Color online) The critical fraction 0j((;/>, Tp, Ti) as a 
function of Ti and 4> for Tp = 2.0 x 10"^. 



V. CRITICAL SCALINGS OF THE VALUE OF 
PLATEAU AND THE PRESSURE 



minimum by using the conjugate gradient technique [39|. 
Then, if the pressure p at the potential energy minimum 
is higher than a threshold value Pth, we increase the vol- 
ume per particle v by Aw. Here, in order to increase 
the volume w, we change the system size L and the po- 
sition Ti of the z-th particle as L x {(i; -I- Az;)/w}^/'^ and 
rj X {(w + At;)/w}^/'^, respectively. After the change of 
the volume, the system is relaxed to the nearest potential 
energy minimum. We repeat the decrease of the volume 
fraction or expand the volume, and relax the system to 
a steady state. Finally, the critical fraction i^j is de- 
termined from the volume per particle wj where the the 
pressure becomes lower than Pth as = Va,v/vj with 
Uav = TT X^i^i '^i I (6./V) . If the pressure at the initial min- 
imum of the potential energy is lower than Pth, we de- 
crease the volume v by Au, relax the system, and repeat 
the decrease and the relaxation until the pressure exceeds 
Pth- Then, we can determine the critical fraction from 
the volume per particle wj where the the pressure exceeds 
Pth. We use Pth = 10"'^ and Aw = 0.0005. We have 
checked that the critical fraction does not change if we 
use Pth = 10"^ and Aw = 0.00005. It is also noted that 
the method to determine the critical fraction is almost 
identical to that in the previous works [1, [^, . 



In Fig. [SI we display the critical fraction (/)j(0, Tp, Tj) 
as a function of (/) and Tp for T\ = 2.0 x lO"-^. Figure [7] 
exhibits 4'j{(t>, Tp, 7i) as a function of the initial temper- 
ature Ti and (f> for Tp = 1.0 x 10"'*. These figures reveal 
that the critical fraction cpj depends on the initial equilib- 
rium state, i. e. Tj and 0, and the quenched state at Tp. 
We note that the existence of the initial state dependence 
has already numerically demonstrated in Ref. [36j . but 
the dependence on the quenched state at Tp within our 
knowledge has not been discussed in any other papers. 



In this section, let us develop the scaling analysis to 
characterize the behavior of MSD and the pressure. This 
section consists of three parts. In the first part, we sum- 
marize some asymptotic relations in the scaling functions. 
In the second part, we briefly introduce the method to 
evaluate the scaling exponents. In the last part, we dis- 
cuss the values of the critical exponents. 

Through our simulation, we have confirmed that that 
the value of plateau mp((/), Tp, Tj) and the pressure 
Tp, Ti) satisfy the scaling laws with the protocol de- 
pendent critical fraction (/)j(0, Tp, 0): 



mp(0,Tp,Ti) = T^-M 



p(^,Tp,Ti) = Tp'-P 



0j(0,Tp,Ti) 
/)j((?!),Tp,Ti) 



(4) 



(5) 



where am, cLp, and b are the critical exponents. Here, we 
assume that the scaling functions M{x) and P(a;) satisfy 



lim M(a;) cx 

a;— >-oo 

lim M{x) cx 

x—>- — oo 

lim P(a;) cx x'^"^'', 
lim P(x) cx 

because of the relations 



lim rrir. 



Fi{ 



lim p = TyF2{4> ~ 4>j), 

Tf-!-0 



for 4> < <p], and 



lim nip = TpP3(0-0j 

Tf^O 

lim p = F4^{4> - 

Tf-!-0 



(6) 
(7) 
(8) 
(9) 

(10) 
(11) 

(12) 
(13) 
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for (p > (pj, where Fi, F2, F3, and F4 are functions 
depending only on — The corresponding scahng 
forms have aheady discussed in terms of the rephca the- 
ory [H, 111]. The similar critical scaling laws are also 
found for the iamming transition for sheared frictionless 
particles ^^U^. 

Figures [5] and IH] show the scaling plots based on Eqs. 
dU and (0) for Tj = 2.0 x lO'^ and 5.0 x IQ-^, respec- 
tively. These figures confirm the validity of Eqs. ^ and 
([5]). Here, we numerically estimate 



10^ 



10" 



10" 



= 0.62, 
= 0.63, 
= 0.64, 
= 0.65, n 
-- 0.66, 5 
= 0.67, 
= 0.68, 
= 0.69, 



I = 0.05 
t = 0.05 
t = 0.05 
5 = 0.05 
t = 0.05 
5 = 0.05 

r = 0.05 

7^ = 0.05 




0.722±0.004, 



0.506±0.004, 



b = 0.471±0. 
(14) 

for different initial temperatures Tj — 2.0 x 10 ^,5.0 x 
10-3, 1.0 X 10-2, 2.0 X 10-2, 3.0 X 10-2, 4.0 x IQ-^, and 
4.0 x 10^2 by using the Levenberg-Marquardt algorithm 
[3^ , where we expand the functional forms of the scaling 
functions as 



004 



ELo^"log(x)" (a;>l), 
ELoBnlogixr ix<l), 

ELo^nlog(a;)" (a:<l). 



(15) 
(16) 



A 



Cn 



and 

the values of the fit- 

A2,A3,A4,A5) 



with fitting parameters 
Dn- Here, we estimate 
ting parameters as {Ao,Ai 
(1.3,-5.2,4.0,-6.2,2.6,-0.8), 

{Bo, Bi,B2, 33,64,35) = 
(2.2, 5.2, 16, -74, -167, -113), (Co, ^1,^2,^3, C4, C5) = 
(-0.9,-8.0,-12,19,-11,4.1), 
{Do,Di,D2,D3,D4,D5) 

(-1.2,0.4,-32,-42,-36,-5.8). This method has 
been used to estimate the critical exponents for the 
iamming transition for sheared frictionless particles 
]40| and for sheared frictional grains [23 |. As shown in 
Figs. [5]and|ni rup and p for different initial temperature 
Ti satisfy the critical scalings with the same critical 
exponents. This indicates that the critical exponents 
are independent of the protocol although the critical 
fraction depends on it as demonstrated in the previous 
section. 

Now, let us estimate the critical exponents in Eqs. (UJ 
and (O by using the previous phenomenological results 
on the jammed soft particles without temperature [1, 
and the unjammed hard spheres [33,|4l|. From Eqs. (jj) 
-(HI, we readily obtain 



lim nip cx |0 — 



for 



< 



lim p cx Tp\(j) — (j>j 

Tp-i-O 



and 

lim m„ cx Tvldi ~ 
Tf^O ' 

lim p cx \(f> — (j)] 

Tp — ^0 



(ap-l)/& 



(a„-l)/b 



(17) 
(18) 

(19) 
(20) 



FIG. 8: (Color online) Scaling plots of mp characterized by 
Eq. © for Ti = 2.0 x IQ-^ and 5.0 x IQ-^ with U = 10^ 
and cj> = 0.62,0.63,0.64,0.65,0.66,0.67,0.68, and 0.69. The 
solid line is the scaling function given by Eq. psp with the 
exponents estimated as Eq. (1141) . The values of the other 
fitting parameters are shown in the text. 




1^^ - h\ I V 

FIG. 9: (Color online) Scaling plots of p characterized by Eq. 
® for Ti = 2.0 x 10"^ and 5.0 x IQ-^ with U = 10^ and p = 
0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, and 0.69. The solid lines 
are the scaling function given by Eq. (|15p with the exponents 
Eq. p4p . The values of the other fitting parameters are shown 
in the text. 



For (/)<(/) J, the pressure may satisfy 
p cx Tf|0- (/)jr^ 



(21) 



as suggested by the free volume theory for hard sphere 
liquids [41]. From the comparison of this equation with 
Eq. ([T^. we obtain 



(22) 



For (/)>(/) J , the pressure might be given by d, Q 



for (b > 



p cx ]0- 0j]. 

From Eqs. (HOI) and we obtain 
Op ^ ^ 



(23) 



(24) 
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For (/><(/) J, MSD for hard sphere hquids is expected 
to satisfy |34] 

mp ~ \(p-h?'^- (25) 
Thus, from Eqs. ^\ and ([25]), we obtain 



(26) 



From Eqs. (gH), (HH), and (gH), we obtain the critical 
exponents 



1 



(27) 



which are not far from the numerical estimated exponents 
presented in Eq. p^ . 



identical critical scaling if we assume that there exists 
only one characteristic length scale, but the values of the 
critical exponents 7 = = 1/2 from the replica theory 
differ from those for nip in Eq. p^ . 

In this paper, we only consider the system with 
Hookean soft-core repulsion given by Eq. ([1]). For a 
system with the interaction potential 



y{n]) = (1 - r^j/a^j) 



A+1/ 



(c^^. 







(29) 



with a exponent A, the scaling of the 
Eq. (1^ is expected to be changed as 

which leads to 



)ressure given by 
li 

(30) 



VI. DISCUSSION AND CONCLUSION 

Now, let us discuss and conclude our results. First, 
we discuss the relationship between our approach and 
the papers by Berthier and Witten [H, [sO]- Second, we 
compare our results with the prediction by the replica 
theory (ssj . Third, we comment on the possibility to 
extend our model to another model of contact. In final, 
we summarize our results. 

The previous papers [1^ H^l demonstrate that the 
structural relaxation time satisfies a scaling relation 
around (/)g = 0.635. We could also reproduce their scal- 
ing in our simulation, though the results are not reported 
in this paper. The scaling relation means that the time to 
escape from the plateau of (f^(i)) for hard sphere liquids 
diverges at (^q ^ but the value of plateau does not exhibit 
any criticality around (pQ because the particles can move 
in the cage even at 4>q. Moreover, the pressure continu- 
ously changes around (pQ because the divergence of the 
relaxation time is not related to the pressure. Hence, 0g 
does not appear in the scaling relations ^ and ([5]). From 
Eq. (|4]), nip for hard sphere liquids becomes zero at 
which indicates that the dynamics of the particles in the 
cage is frozen at j33 |. 

In Ref. [33|, the replica analysis is used for the ex- 
planation of the jamming transition on the temperature- 
density plane for harmonic spheres. They derived the 
identical scaling relation ([5]) for the pressure with the 
critical exponents = 6 = 1/2 corresponding to Eq. 
([27)1 and our numerical results in Eq. (fT4l) . In addition, 
they suggested a critical relation 



b 



(31) 



A* = T'^A' 



(28) 



for the optimal cage size A* with critical exponents 
7 = 1/ = 1/2. Because both of A* and nip are the charac- 
teristic length, it may be reasonable that they satisfy the 



From Eqs. (j^ . (j^ . and (PT|) . the critical exponents 
for the system with the potential given by Eq. (|29|) are 
expected to be 



A 1 , , 

a„ = ^^, b=^^. (32) 



2A + 2' 



A + 



A + 1 



The similar dependence of the critical exponents is con- 
firmed in the sheared granular systems |18| . 

In conclusion, we have numerically investigated critical 
behavior of VP of MSD and the pressure for soft repul- 
sive particles after quench near the jamming transition 
point. We verify the existence of the critical scaling re- 
lations Q and ([5]), and numerically evaluate the critical 
exponents and the critical fraction. The critical fraction 
exhibits the protocol dependence, while the critical ex- 
ponents are independent of the protocol, which are close 
to the estimation Eq. (P?]) in terms of the combination 
of the existing arguments. 
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